# This file produces tables 3 and 4
# Set working directory
setwd("/Users/ericguntermann/Documents/Papers/Representation of Party Preferences/Replication")

set.seed(123)
library(foreign)

## Criterion 1 (proportional)
criteria <- read.dta("criteria.dta")
criteria <- criteria[-which(criteria$countryyear %in% c("Thailand 2001", "Thailand 2011", "Japan 1996", "Japan 2013")),]
d1 <- with(criteria, list(propfic=propfic, proportional=proportional, gdppercap=gdppercap, freedomhouse=freedomhouse))
d1$N <- length(unique(criteria$countryyear))

m1 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.proportional*proportional[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.proportional ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

p1 <- c("b", "b.proportional", "b.gdppercap", "b.freedomhouse")

i1 <- function(){
  list("b"=rnorm(1), "b.proportional"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}


## Criterion 1 (gallagher)
library(foreign)
dat <- read.dta("criteria.dta")
d2 <- with(dat, list(propfic=propfic, gallagher=log(gallagher), gdppercap=gdppercap, freedomhouse=freedomhouse))
d2$N <- length(unique(dat$countryyear))
m2 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.gallagher*gallagher[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.gallagher ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

i2 <- function(){
	list("b"=rnorm(1), "b.gallagher"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}
p2 <- c("b", "b.gallagher", "b.gdppercap", "b.freedomhouse")


## Criterion 1 (mdm)
library(foreign)
criteria <- read.dta("criteria.dta")
d3 <- with(criteria, list(propfic=propfic, mdm=log(mdm), gdppercap=gdppercap, freedomhouse=freedomhouse))
d3$N <- length(unique(criteria$countryyear))

m3 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.mdm*mdm[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.mdm ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p3 <- c("b", "b.mdm", "b.gdppercap", "b.freedomhouse")


i3 <- function(){
  list("b"=rnorm(1), "b.mdm"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}

## Criterion 1 (proportional) with number of parties
criteria <- read.dta("criteria.dta")
criteria <- criteria[-which(criteria$countryyear %in% c("Thailand 2001", "Thailand 2011", "Japan 1996", "Japan 2013")),]
d4 <- with(criteria, list(propfic=propfic, proportional=proportional, gdppercap=gdppercap, freedomhouse=freedomhouse, ngov=ngov))
d4$N <- length(unique(criteria$countryyear))

m4 <- function(){
  for(i in 1:N){
    propfic[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.proportional*proportional[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]+ b.ngov*ngov[i]
  }
  b ~ dnorm(0,0.01)
  b.proportional ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  b.ngov ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

p4 <- c("b", "b.proportional", "b.gdppercap", "b.freedomhouse", "b.ngov")

i4 <- function(){
  list("b"=rnorm(1), "b.proportional"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1), "b.ngov"=rnorm(1))
}


d <- list(d1,d2,d3,d4)
p <- list(p1,p2,p3,p4)
m <- list(m1,m2,m3,m4)
ini <- list(i1,i2,i3,i4)


library(foreach)
library(doParallel)
library(R2jags)
registerDoParallel(4) 
cl <- makeCluster(4)
registerDoParallel(cl)
results <- foreach(i=1:4, .packages = c('R2jags')) %dopar% jags(data=d[[i]], model=m[[i]], parameters.to.save=p[[i]], inits=ini[[i]], n.iter=100000, n.burnin=50000, n.thin=10, n.chains=2)

results1 <- as.mcmc(results[[1]])
results2 <- as.mcmc(results[[2]])
results3 <- as.mcmc(results[[3]])
results4 <- as.mcmc(results[[4]])


# Convergence diagnostics

library(ggmcmc)
p1 <- ggs_geweke(ggs(results1)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 1")
p2 <- ggs_geweke(ggs(results2)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 2")
p3 <- ggs_geweke(ggs(results3)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 3")
p4 <- ggs_geweke(ggs(results4)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 4")
library(gridExtra)

pdf("criterion1_conv.pdf")
grid.arrange(p1,p2,p3,p4)
dev.off()

## GET OUTPUT
#Model 1 (Proportional)
coefs1 <- as.data.frame(as.matrix(results1))
b <- coefs1[,"b"]
b.proportional <- coefs1[,"b.proportional"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co1 <- c(mean(b), mean(b.proportional), NA, NA, mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b), sd(b.proportional), NA, NA, sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.proportional)>0, mean(b.proportional>0),mean(b.proportional<0)), NA, NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 2 (Gallagher)
coefs1 <- as.data.frame(as.matrix(results2))
b <- coefs1[,"b"]
b.gallagher <- coefs1[,"b.gallagher"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]


co2 <- c(mean(b), NA, mean(b.gallagher), NA, mean(b.freedomhouse), mean(b.gdppercap))
sd2 <- c(sd(b), NA, sd(b.gallagher), NA, sd(b.freedomhouse), sd(b.gdppercap))
p2 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, ifelse(mean(b.gallagher)>0, mean(b.gallagher>0),mean(b.gallagher<0)), NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 3 (MDM)
coefs1 <- as.data.frame(as.matrix(results3))
b <- coefs1[,"b"]
b.mdm <- coefs1[,"b.mdm"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co3 <- c(mean(b), NA, NA, mean(b.mdm), mean(b.freedomhouse), mean(b.gdppercap))
sd3 <- c(sd(b), NA, NA, sd(b.mdm), sd(b.freedomhouse), sd(b.gdppercap))
p3 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, NA, ifelse(mean(b.mdm)>0, mean(b.mdm>0),mean(b.mdm<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))

# Table3
library(xtable)
model.effects <- data.frame(co1,sd1,p1,co2,sd2,p2,co3,sd3,p3)                                                                                                                
rownames(model.effects) <- c("Intercept", "Proportional", "Gallagher*", "MDM*", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P", "Mean", "SD", "P","Mean", "SD","P")
xtab <- xtable(model.effects, caption="Criterion 1", digits=c(2))
sink("criterion1.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)

#Model 4 (Proportional with ngov)
coefs1 <- as.data.frame(as.matrix(results4))
b <- coefs1[,"b"]
b.proportional <- coefs1[,"b.proportional"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]
b.ngov <- coefs1[,"b.ngov"]

co1 <- c(mean(b), mean(b.proportional), mean(b.ngov), mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b), sd(b.proportional), sd(b.ngov), sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.proportional)>0, mean(b.proportional>0),mean(b.proportional<0)), ifelse(mean(b.ngov)>0, mean(b.ngov>0), mean(b.ngov<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))

# Table4
library(xtable)
model.effects <- data.frame(co1,sd1,p1)                                                                                                                
rownames(model.effects) <- c("Intercept", "Proportional", "Number of parties", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P")
xtab <- xtable(model.effects, caption="Criterion 1", digits=c(2))
sink("criterion1b.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)

